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Abstract 

A C++ code for the solution of the relativistic Hartree-Bogoliubov theory in coordinate 
space is presented. The theory describes a nucleus as a relativistic system of baryons and 
mesons. The RHB model is applied in the self-consistent mean-field approximation to the 
description of ground state properties of spherical nuclei. Finite range interactions are 
included to describe pairing correlations and the coupling to particle continuum states. 
Finite element methods are used in the coordinate space discretization of the coupled 
system of Dirac-Hartree-Bogoliubov integro-differential eigenvalue equations, and Klein- 
Gordon equations for the meson fields. The bisection method is used in the solution of the 
resulting generalized algebraic eigenvalue problem, and the biconjugate gradient method 
for the systems of linear and nonlinear algebraic equations, respectively. 
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PROGRAM SUMMARY 


Title of program: spnRHBfem.cc 

Catalogue number: . 

Program obtainable from: 

Computer for which the program is designed and others on whieh it has been tested: any 
Unix work-station. 

Operating system: Unix 

Programming language used: C-|—|- 

No. of lines in eombined program and test deek: 

Keywords: relativistic Hartree-Bogoliubov theory, mean-field approximation, spher¬ 
ical nuclei, pairing, Dirac-Hartree-Bogoliubov equations, Klein-Gordon equation. Fi¬ 
nite Element Method, bisection method, classes 


Nature of physical problem 

The ground-state of a spherical nucleus is described in the framework of relativistic 
Hartree-Bogoliubov theory in coordinate space. The model describes a nucleus as a 
relativistic system of baryons and mesons. Nucleons interact in a relativistic covari¬ 
ant manner through the exchange of virtual mesons: the isoscalar scalar cj-meson, 
the isoscalar vector u;-meson and the isovector vector p-meson. The model is based 
on the one boson exchange description of the nucleon-nucleon interaction. Pairing 
correlations are described by finite range Gogny forces. 

Method of solution 

An atomic nucleus is described by a coupled system of partial integro-differential 
equations for the nucleons (Dirac-Hartree-Bogoliubov equations), and differential 
equations for the meson and photon fields (Klein-Gordon equations). A method 
is presented which allows a simple, self-consistent solution based on finite element 
analysis. Using a formulation based on weighted residuals, the coupled system of 
Dirac-Hartree-Bogoliubov and Klein-Gordon equations is transformed into a gener¬ 
alized algebraic eigenvalue problem, and systems of linear and nonlinear algebraic 
equations, respectively. Finite elements of arbitrary order are used on adaptive 
non-uniform radial mesh. The generalized eigenvalue problem is solved in narrow 
windows of the eigenparameter using a highly efficient bisection method for band 
matrices. A biconjugate gradient method is used for the solution of systems of linear 
and nonlinear algebraic equations. 




Restrictions on the complexity of the problem 

In the present version of the code we only consider nuclear systems with spherical 
symmetry. 


LONG WRITE-UP 


1 Introduction 

Relativistic mean-field theory has been extensively applied in calculations of nuclear mat¬ 
ter and properties of finite nuclei throughout the periodic table. The theory provides 
a framework for describing the nuclear many-body problem as a relativistic system of 
baryons and mesons [Q, In the self-consistent mean-field approximation, detailed 

calculations have been performed for a variety of nuclear stucture phenomena (for a recent 
review see |^). More recently, the relativistic mean-field model has been also applied in 
the description of the structure of exotic nuclei with extreme isospin values. This new field 
includes many interesting phenomena: extremely weak binding of the outermost nucleons, 
coupling between bound states and the particle continuum, regions of neutron halos with 
very diffuse neutron densities, large spatial dimensions and the existence of the neutron 
skin. Major modifications of shell structures in drip-line nuclei have been predicted, as 
well as modifications in the onset and evolution of collectivity. For most phenomena 
along the line of stability, non-relativistic models and the relativistic framework predict 
very similar results, although relativistic mean-field models provide a more economical 
description. For drip-line nuclei, on the other hand, one also expects differences in the 
predictions of non-relativistic and relativistic models, especially in the treatment of the 
spin-orbit interaction |Q]. 

In drip-line nuclei the Fermi level is found close to the particle continuum. The lowest 
particle-hole or particle-particle modes are often embedded in the continuum, and the 
coupling between bound and continuum states has to be taken into account explicitly. 
In the mean-field approximation the most important residual interaction is the pairing 
force. Therefore, in the description of ground-state properties, it is essential to include 
mean-field and pairing correlations simultaneously. The Relativistic Hartree-Bogoliubov 
(RHB) theory in coordinate space, which is an extension of non-relativistic HFB-theory 
provides such a unified description. In particular, it includes the scattering of nucleonic 
pairs from bound states to the positive energy continuum. The RHB theory has recently 
been applied in the description of ground-state properties of Sn and Pb isotopes |^, using 
an expansion in a large oscillator basis for the solution of the Dirac-Hartree-Bogoliubov 
equations. In many applications an expansion of the wave functions in an appropriate 
oscillator basis of spherical or axial symmetry provides a satisfactory level of accuracy. 
However, in the case of drip-line nuclei, the expansion in the localized oscillator basis 
presents only a poor approximation to the continuum states, and the convergence of such 
expansions is too slow and is not uniform. Examples are exotic phenomena such as neutron 
halos and neutron skins. In order to correctly describe the coupling between bound and 
continuum states, the Dirac-Hartree-Bogoliubov equations have to be solved in coordinate 
space. Discretization in coordinate space provides also the advantage that exotic shapes 
and/or large deformations can be described, without preparing a basis specihc to each 
deformation. Recently, a fully self-consistent RHB model in coordinate space has been 
used to describe the two-neutron halo in ^^Li Q. However, only a density dependent force 


of zero range has been used in the pairing channel. In general, it is assumed that a finite 
range interaction would provide a more realistic description of pairing correlations. In 
the present article we describe a C++ code which can be used to calculate ground state 
properties of spherical nuclei, within the framework of RHB theory in coordinate space, 
with finite range forces in the pairing channel. 

A convenient procedure for the coordinate space discretization of the Dirac-Hartree- 
Bogoliubov equations is provided by Finite Element Methods (FEM) |^, ^, 1+, n|. In 
Refs. [1^, 1R| we have applied EEM to the solution of the coupled system of relativistic 
mean-field equations in the description of a one-dimensional slab of nuclear matter and 
of spherical doubly closed-shell nuclei. We have investigated the applicability of FEM 
in the calculation of bound and continuum eigenstates of the Dirac equation. Since the 
spectrum of the hamiltonian of the hyperbolic Dirac equation is not bounded from below, 
finite element methods cannot be applied in the variational formulation. The method of 
weighted residuals produces element matrix integral definitions that would be identical 
to those obtained from a variational form, if one existed. Our analysis has shown that 
FEM provide very accurate solutions for the relativistic eigenvalue problem in the self- 
consistent mean-field approximation. In the present work we extend the model to include 
pairing correlations in the framework of RHB theory. 

The article is organized as follows. In Sec. 2 the Relativistic Hartree-Bogoliubov theory 
is described. Dirac-Hartree-Bogoliubov equations for a system with spherical symmetry 
are derived. The finite element analysis is described in Sec. 3. In Sec. 4 we present 
some illustrative calculations and discuss the quality of our approximations and numerical 
results. The structure of the C++ code is described in Sec. 5. 


2 The relativistic Hartree-Bogoliubov equations 


The Hartree-Eock-Bogoliubov (HFB) theory provides a unified description of mean-field 
and pairing correlations in nuclei |M|. Independent quasiparticles are introduced and the 
ground state of a nucleus |<I> > is represented as the vacuum with respect to these quasi¬ 
particles. The quasi-particle operators are defined by a unitary Bogoliubov transformation 
of the single-nucleon creation and annihilation operators. The generalized single-particle 
hamiltonian of HFB theory contains two average potentials: the self-consistent field F 
which encloses all the long range ph correlations, and a pairing field A which sums up 
the pp-correlations. The expectation value of the nuclear hamiltonian < <I>|.H|<I> > can be 
expressed as a function of the hermitian density matrix p, and the antisymmetric pairing 
tensor k. The variation of the energy functional with respect to p and k produces the 
single quasi-particle Hartree-Fock-Bogoliubov equations (for details of the derivation we 
refer to [1R 


(h-X A \(Uk\ 
l-A* -h + x)\vj 



( 1 ) 


HFB-theory, being a variational approximation, results in a violation of basic symmetries 
of the nuclear system, among which the most important is the nonconservation of the 
number of particles. In order that the expectation value of the particle number operator 
in the ground state equals the number of nucleons, equations contain a chemical 
potential A which has to be determined by the particle number subsidiary condition. The 
column vectors denote the quasi-particle wave functions, and Ek are the quasi-particle 
energies. 









The relativistic extension of the HFB theory is descibed in Ref. |T^ . In the Hartree ap¬ 
proximation for the self-consistent mean field, the Relativistic Hartree-Bogoliubov (RHB) 
eqnations read 


(hD-m-X A \ (Uk\ ^ p (Uk\ 

V -A* -hD + m + x)\Vk) ^\Vk)- 


( 2 ) 


where h^j \s the single-nncleon Dirac hamiltonian and m is the nucleon mass. The 
RHB equations are non-linear integro-differential equations. They have to be solved self- 
consistently, with potentials determined in the mean-field approximation from solutions 
of Klein-Gordon equations for mesons |13] 


[-A + ml]a{r) 

= -9a Ps (r) - 92 (r) - 53 <7^ (r) 

(3) 

-A + ml] a;°(r) 

= -9ujPvir) 

(4) 

[-A + ml] p°(r) 

= -9pP3{r) 

(5) 

- AA°(r) 

= epp(r). 

(6) 


for the sigma meson, omega meson, rho meson and photon field, respectively. The spatial 
componenets cn, p, and A vanish due to time reversal symmetry. Because of charge 
conservation, only the 3-component of the isovector rho meson contributes. The source 
terms in equations (|3|) to @ are sums of bilinear products of baryon amplitudes 


Ps 

= E uVu, 

Ek>0 

(7) 

Pv 

= E 

Ek>0 

(8) 

P3 

= E ^kr3Vk, 

Ek>0 

(9) 

Pern 

Ek>0 

(10) 


where the sums run over all positive energy states. For M degrees of freedom, for example 
number of nodes on a radial mesh, the HB equations are 2M-dimensional and have 2M 
eigenvalues and eigenvectors. To each eigenvector {Uk,Vk) with eigenvalue E^, there 
corresponds an eigenvector with eigenvalue —E^. Since baryon quasi-particle 

operators satisfy fermion commutation relations, it is forbidden to occupy the levels E^ 
and —Ek simultaneously. Usually one chooses the M positive eigenvalues Ek for the 
solution that corresponds to a ground state of a nucleus with even particle number. 

The system of equations @, and to is solved self-consistently in coordinate 
space by discretization on the finite element mesh. In the coordinate space representation 
of the pairing field A in (^), the kernel of the integral operator is 

Aabir, r') = -^ X! ^abcdir, r')Kcd(r, r'). (11) 

^ c,d 

where a, b, c, d denote all quantum numbers, apart from the coordinate r, that specify 
the single-nucleon states. Uabcd(r, rO matrix elements of a general two-body pairing 
interaction, and the pairing tensor is defined as 

Kcd(r,r'):= C/4(r)V'rffc(r')- 

Ek>0 


c: 


( 12 ) 




The integral operator A acts on the wave function 14 (r): 

(A14)(r) = ^ f (fr'Aab{r,r')Vbk{r'). (13) 

b 

The eigensolutions of Eq. (2) form a set of orthogonal (normalized) single quasi¬ 
particle states. The corresponding eigenvalues are the single quasi-particle energies. The 
Bogoliubov transformation from the single-particle coordinate basis of (5-functions to the 
basis of quasi-particle states is given by the matrix Wkr ■= (^)! ^ k are 

column and row indices, respectively. In the self-consistent iteration procedure we work in 
the basis of quasi-particle states. The self-consistent quasi-particle eigenspectrum is then 
transformed into the canonical basis of single-particle states. The canonical basis is defined 
to be the one in which the matrix Rkk' = (14(r)|14/(r)) is diagonal. The transformation to 
the canonical basis determines the energies and occupation probabilities of single-particle 
states, that correspond to the self-consistent solution for the ground state of a nucleus. In 
order to determine the canonical basis, we have two possibilities: either diagonalize the 
density matrix, or diagonalize the matrix 

Rkk' ■■= mr)\Vk>{r)) (14) 


Although both methods are in principle equivalent, for numerical reasons we chose the 
second method, i.e. we diagonalize the matrix Rkk'- Because of the truncation in quasi¬ 
particle space, the dimension of the matrix Rkk’^ which is a matrix in quasiparticle space, 
is considerably smaller than the dimension of the density matrix in coordinate space. 

The transformations from the general single-particle basis to the basis of quasi-particle 
states, and the canonical basis are illustrated in Fig. I. The quasi-particle operators and 
basis states are defined 


= / (Rr 


U*k,r 


Uk,rJ U^(r) 


The operators in the canonical basis are 


= / err 


V’fc(r) 0 W 

0 rk{r)J\^Hr) 



and the transformation from the quasi-particle to the canonical basis reads 


( \ ( '^kk> Vkk’ \ f ak'\ 

V 4 / V ~'^kk’ Ukk' ) V / ’ 


V'fc(r) = '^Ck'k^v,k'i'^) RR 

k' 


The matrix representation of the unitary Bogoliubov transformation W can be decom¬ 
posed into a product of three matrices 0 


W = 


D 0 

0 D* 


u Y 
-V u 


C 0 

0 c* 


(18) 


The diagonalization of the matrix Rkk' '■= (I4(r)|I4'(r)) produces the unitary matrix C. 
Columns of C are eigenvectors of Rkk'- The matrices U and V are constructed from the 


eigenvalues vf, of Rkk’ 
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o 
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Uk 
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Vk 
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where = yl — and N is the number of solutions. The matrix C transforms the 
basis of quasi-particle states into the canonical basis 

V'fe(r) ='^Ck'kVk'{r) 
k' 

The matrix D := [^^(r)] represents a transformation between two single-particle bases. 
E := diag(Efc) defines the diagonal matrix of quasi-particle energies. In the basis of 
quasi-particle states the hamiltonian matrix has the diagonal form diag(E,—E). The 
hamiltonian matrix in the single-particle canonical basis is defined by the transformation 

VWC OWE 0 wet OWtJ -V\ 

■ \-Y \j)\0 CV^vO -eJI^O lyVt Ut j • 

The single-particle energies correspond to the diagonal matrix elements = Hnn + A, 
where A denotes the chemical potential. In Fig. 2 we display a schematic single-nucleon 
spectrum in the relativistic mean-field potential of a finite nucleus. On the left hand 
side the eigenspectrum of a Dirac hamiltonian is shown. The single-particle hamiltonian 
corresponds to the average mean-field potential, and the Dirac equation is solved in the 
Hartree mean-field, and no-sea approximations. A Dirac gap is observed between states 
of negative and positive energy. In vacuum, this gap equals two times the nucleon mass. 
In a nucleus the Dirac gap extends between the sum and the difference of the scalar 
sigma-meson potential and vector omega-meson potential. The sum and the difference 
are given relative to -\-m and —m, respectively. In the center of Fig. 2 the eigenspectrum 
of the Dirac hamiltonian is shifted by the nucleon mass m. As a result, bound single¬ 
nucleon levels have negative energies, while the positive energy domain contains only 
single-nucleon continuum states. The single quasi-particle spectrum which results as a 
solution of the relativistic Hartree-Bogoliubov equations is shown on the right hand side 
of Fig. 2. The number of solutions is two times the number of physical states. As already 
described, for each eigenvector {Uk, 14 ) with energy Ek, the corresponding state {V^, U^) 
is found at —Ek- 

In practical calculations the Dirac-Hartree-Bogoliubov and Klein-Gordon equations 
are discretized on a finite domain V in coordinate space (indicated by Vmax in Fig. 2), 
and the generalized eigenvalue problem is solved in the window £ := [0,Emax] of the 
eigenparameter E. The domain E 0 £ is indicated by the shaded area in the right hand 
side spectrum. By increasing the coordinate space domain the HB spectrum becomes 
denser, and thus provides a better approximation for the continuum. Larger values of 
Emax take into account couplings to highly excited quasi-particle states. V and £ should 
be chosen in such a way that the resulting densities do not depend on their precise values. 
In particular, E^ax has to be larger than the absolute value of the depth of the potential 
well. 



( 20 ) 



In the present version of the code we only consider single closed-shell nuclei, i.e. sys¬ 
tems with spherical symmetry. The fields (7{r), uj^(r), p^{r), and A^{r) depend only on 
the radial coordinate r. The nucleon spinors Uk (14) in (|^ are characterized by the an¬ 
gular momentum j, its 2 ;-projection m, parity vr and the isospin ta = for neutron and 
proton. We combine the two Dirac spinors ?7fc(r) and 14(r) to form a super-spinor 


^fc(r) := 


Wkir)) 


( 22 ) 


where 




(23) 


g{r) and /(r) are radial amplitudes, Xt is the isospin function, the orbital angular mo¬ 
menta I and I are determined by j and the parity vr 


and 


I = 


I = 


j -|- 1/2 for vr = (—1)'^+^/^ 

j — 1/2 for vr = ( — 1)-^“^/^ ’ 


j — 1/2 for vr = (— 

j 1/2 for vr = 


Hjim is the tensor product of the orbital and spin functions 


(24) 


(25) 




rris^mi 


(26) 


It will be useful to define a single angular quantum number k as the eigenvalue of the 
operator (1 -|- d • 1) 

(1 T d • i ) Hkj-tti — ki ^K,mi (^^) 

K = ±{j + l/2) for j = 1^1/2. (28) 

K = ±1, ±2, ±3,..., and the Dirac HB equations are solved in coordinate space for each 
value of K. The equations for the radial amplitudes gu{v){^) fu{v){''") derived 

from Eq. (|^. The radial single-particle Dirac Hamiltonian reads (see also appendix A) 

h^{r) -m:= + r~^) - <7iKr~^ -h (as - l 2 )m -b cr^Sir) -\- l 2 Vb(r), (29) 


where the scalar and vector potentials are 

S{r)=g„a{r), and H°(r) = gujUJ°{r) +gpT 3 pl{r) + ^^^^ A°(r). (30) 

iTj (i=l,2,3) are the Pauli matrices. The integral operator of the pairing interaction takes 


the form 

OO 

A(r) = J dr'r''^A{r,r') 

0 

The kernel of the integral operator is defined 

^i^ir,r') = ^^{ra,ra'\V\rd,r'd')j^Ka,a'{r,r') 

d,d' 


(31) 


(32) 


O 



where r and r' denote radial coordinates, o, a', a and a' are quantum numbers that com¬ 
pletely specify single-nucleon states: (n, l,j, m) or (n, k, m), J and M are the total angular 
momentum of the pair, and its z-projection, respectively. For the pairing interaction we 
use the finite range Gogny force 

(ri -^2)^ 

VPP(1,2}= e {Wi +BiP^ - HiP^ - MiP^P^). (33) 

i=l,2 

We only consider contributions from J = 0 pairs to the pairing matrix elements. The 
kernel of the integral operator can then be written 

= r')Kf,{r, r). (34) 


Details on the calculation of two-body matrix elements of the Gogny force are given in 
Appendix B. The pairing tensor is calculated 


K, 


:Xr,r') = Y.2\ 


K\ 


sIa {r') 0 

0 


(35) 


where, for a quantum number k, the sum runs over all solutions in the specified energy 
interval 0 < < T'max- If we define ‘h^(y)(r) := ( 5 t/(y)(r),/[/(v')(r))^, the radial Dirac- 

Hartree-Bogoliubov equations read 


poo 

{hoir) — m — + / (irV^A(r,/)4>y(/) = i?<l>[/(r) 

Jo 

poo 

{—hDir)+m + X)^v{r)+ / (irV^A(r, r')$t/(r') = i?<l>y (r) 

Jo 

The meson and photon fields are solution of the Klein-Gordon equations 


(36) 


, . 2^ 1(1 +1) 9 . . , 

(dr ^dr+ ^2 +ml)a(r) 

= -9a Ps (r) - 92 <7^ (r) - 53 (r) 

(37) 

. ^dr+ *'^2 

= -9ujPv(r) 

(38) 

( ^dr+ ^^2 

= -9pPo(r) 

(39) 

(-9'-^5, + ^^^)A0(r) 

= epp(r). 

(40) 


where I = 0 for spherical nuclei. In the following we describe the computer code used to 


solve the system of equations (36) to (pOD in a self-consistent iteration scheme. 


3 Finite Element discretization of the radial equa¬ 
tions 

The radial equations (^^ - ( ^0| ) can be written in the general form 

F(A,r,cr(r),a;°(r),p°(r), A°(r), A(r), A)4'i(r) = ei4>i(r) (41) 


n 








Ra{r,a{r))a{r) 

= s<^(4>i(r),.. 

■•,^A(r)) 

(42) 

R^{r)u)^{r) 

= s^(4>i(r),.. 


(43) 

Rp{r)p°{r) 

= Sp($i(r),.. 


(44) 

Rc{r)A^{r) 

= sc(^i(?’),- 


(45) 


where if is a 4 x 4 matrix operator defined as 


/ S{r 

) + V{r) - A 

dr 

- {Ki - l)r~^ 

—dr 

- {Ki + l)r~^ 

—2m 

- S{r) + V{r) - A 


^(99) (r) 


0 

V 
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A(//)(r) 

and 





Ra '■ = 

1 

1 

2 r~^ dr + l{l + 1) 


Ruj '■ = 

1 

1 

2r-^dr + l{l + l) 


Rp '■= 

-d^- 

2r—ldr -I- -I- 1) 


Rc ■■= 

1 

1 

2r~^dr -I- l{l -I- 1) 


H : = 


A(s9)(r) 

A(//)(r) 

0 

-5(r)-F(r) + A 

-dr + {Ki - l)r“^ 

dr + {Ki + l)r~^ 

2m + S{r) — V{r) + X/ 


+ml + g2a{r) + g3a‘^{r) 

(47) 

-2 1 2 

(48) 

-2 1 2 
■ +mp 

(49) 

-2 

(50) 


The source terms on the right hand sides of these equations are defined as ■= —ga Psif), 
Su) '■= guiPvii"), Sp := QpPsir), and sc ■= e/3em(?’)- The method of finite elements is used 
to discretize the system of Dirac-Hartree-Bogoliubov and Klein-Gordon equations on a 
spherical mesh. For the nucleon spinor we define T(r) := <l>y(r)^)^, and use the 

FEM ansatz 

^{r) = J2XpNp{r) (XpeK^), (51) 

p 

where the coefficients Xp are four-component vectors, and Np(r) denotes Lagrange shape 
functions of arbitrary order |j^ . The index p enumerates nodes on a finite element mesh 
for a spherical box (rmin = 0 to rmax)- fp denotes the radial coordinate of the node p. 
The Lagrange shape functions Np satisfy the property Np{rpi) = 5ppi. Therefrom the 
coefficients Xp in the ansatz ( |5l|) correspond to the amplitudes T(rp) of the solution: 
^P = The Dirac-Hartree-Bogoliubov equations 

can be formally written 

b{^) = f{^,x) (52) 

where Z) is a differential operator. For any approximate solution 'I'(x), the residual error 
is defined as 

i?(x) :=D(§)-/(T,x). (53) 

If (1^ has to be solved on a compact domain D C R”, where 'I'(x) = g{x) for xedQ., and 
g{x) is defined on the boundary dO- of D, the method of weighted residuals can be used 
to define the weak form of the boundary value problem 


J R{x) w{x) cT'x = 0 for all weight functions u;(x). (54) 

n 


The method of weighted residuals necessitates that a weighted integral of the residual 
vanishes. In the Galerkin method the choice of weight functions is Wp{x) = Np{x), where 


.( 46 ) 


1 n 





Np{x) are shape functions that dehne the FEM ansatz of the solution. As we have shown 


in Ref. [|l^ , the choice Wp{r) = r^Np{r) produces symmetric stiffness matrices for the 
radial equations. Using the standard representation for the Pauli matrices, the radial 
relativistic Hartree-Bogoliubov equations (^) can be written in matrix form 

{dr + r~^) <J3 ( 8 ) (TsfTi - Kr~^ (T 3 (g) (Ti + m (T 3 (g) ((73 - I 2 ) + S{r) (73 (g) (73 + 


Voir) (73 (g) I 2 + (7i (g) l 2 A(r) - A (73 


'l'(r) — E l 4 'I'(r) = 0 


(55) 


Eq. (^) represents a system of four coupled integro-differential equations. Eor the ma¬ 
trices which define the structure of these equations we use the notation 


Ai := 0-3 (g)o-3cri, 

A 4 := as (g) fJ 3 , 


A 2 := 0-3 (g) fJi, A 3 := 0-3 (g) ( 0-3 - I 2 ), 
A5 := (73 (g) I2, Ag := (7i (g) I2; 


(56) 


The method of weighted residuals transforms the Dirac-Hartree-Bogoliubov equations into 
a hnite system of algebraic equations. The system forms a generalized eigenvalue problem 
AX = eN A with global stiffness matrices 

A = {wp/\{dr + r~^) Ai — Kr~^ A 2 + m A 3 — A A 5 -|- S{r) A 4 -|- 

Uo(r)A5 + A6A(r)|iVp) (57) 

and 

N = ('u;p/|l4|iVp). (58) 

The number of equations in the system is thus four times the number of nodes on the 
finite element mesh. The eigenvalue matrix equation reads 


Si (g) Ai -I- S 2 <g) Ai — KS 2 (g) A 2 + mS 3 (g) A 3 + S 4 (g) A 4 + S 5 (g) A 5 -1- Sg ® Ag — AS 3 (g) Ac 


X 


= e 


x{59) 


where x is a vector with components {x^i’-^\ x[^’-^\ x[^^\ ...,Xn^\xn’-^\xn^\xn^"*)'^. 

The global stiffness matrices Si to Sg correspond to the various operator terms in Eq. 


51 

5 2 

53 

5 4 

55 
Se 


{wp'{r)\dr\Np{r))^, 

{wp'{r)\r-^\Np{r))^, 

{wp'{r)\Np{r))^, 

{wp>{r)\S{r)\Np{r))'\, 

{wp'{r)\V{r)\Np{r))], 

{wp'{r')\A{r',r)r^\Np{r)) 


(60) 


Details of the calculation and construction of the matrices Si to S 5 have been described in 
Appendix A of Ref. |13|. The matrix elements of Sg have to be calculated for the nonlocal 
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pairing interaction A{r',r). Using the same notation as in Ref. []R|, the matrix elements 
of the local stiffness matrix in the reference element representation are 


S6°^= [(iVg'k''A(/,r)r2|iV^) 


/ dp' / dpNq>{p){(p/ + A^p) 


/\3m—1 


A(a(V + ^CpT, a{Cp + ACpDiCp + ACp)^™"^ A,(p) 


(61) 


For the FEM discretization of the Klein-Gordon equations we use the ansatz 


a{r) = ^apNp{r) 

P 

(62) 

u;\r) = Y,u;^pNp{r) 

P 

(63) 

P^ir) = ^plNp{r) 

P 

(64) 

A\r) = Y,AlNp{T). 

(65) 


p 


where the node variables ap, w®, Pp and correspond to field amplitudes at the mesh 
point p. For the Klein-Gordon equations we use the same type of shape functions Np(r), 
and the same mesh as in the FEM discretization of the RHB equation (|55|). In this way 
we obtain the following algebraic equations for the mean field amplitudes of the meson 
fields 


+ l{l + 1 ) + ml 

5- + l{l + 1 ) ^ 2 - + ml S‘^ 


Sl + l{l + l)SP + mlS^ 


Sf + l{l + 1) 52^^° 


a = 


(P = r-< 3 ) 

^0 _ ^em) _ 


The node variables 


/ lO nO 
^'P ? Pp 

and Mp are grouped into the vectors 

a;° = {ioi, p 

>0 _ 

(P?,-, 

pIY, 

10 

= (-4i,.-,-40)^, and 

Si = 


= s{ = 

--st 

= 

{wp,{r)\dl + 2r-^ dr\Np{r)), 

^2 = 

QLO 

= S^2 = 

- ^2 

= 

{wp,{r)\r-lNp{r)), 

^3 = 

QiO 

*^3 

= SP = 

- O 3 

= 

{wp'{r)\Np{r)), 


SI = {wpi{r)\g 2 a{r) + g‘icr{rf Np{r)). 
The components of the right hand side vectors are defined as 

rp' = -9a{wpPr)\p^{r)), 

= gLo{wp'{r)\pPr)), 

^p^ = 9p{wp'{r)\p's,{r)), 

T'pP'^ = e{wp,{r)\pem{r)). 


( 66 ) 

(67) 

( 68 ) 
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(73) 


( 74 ) 


1 o 
















Details on the calculation of the stiffness matrices 8^,82,83, 84,8^, 82,83,8]*,8^ 
8 ^°,8^°, and the right-hand side vectors ,r^'"),^^em)^ given in Appendix B of 


cP 

^ 3 ) 


Ref. 113]. Next we discuss the structure and occupation pattern of the global stiffness 
matrices of the RHB equation and describe the inclusion of boundary conditions. The 
boundary conditions for the spherical symmetric case are 


= 0) = 0 and 
gi^) (r = 0) = 0 and 

g{U) _ Q 

gi'^) (r = 0) = 0 and 


f^^\r = 0) = 0 for K = —1 
gi^)(^r = 0)=0 for K=+l 
fiU)(j. = 0) = 0 for |k| > 1 
= 0) = 0 for |«;| > 1 


(75) 


and 


(?’max) = 0 and = 0 for all K. 


aiU),r 

y V' max 


(76) 


In Figs. 3 (a) and (b) we display the occupation pattern of the global stiffness matrics 
A and N, respectively. The occupation pattern corresponds to first order finite elements 
(linear shape functions). Shaded squares indicate the positions of non-zero matrix ele¬ 
ments. The occupation pattern displays a dominant block diagonal band structure. This 


structure results from the contributions of the local terms in (55). Each diagonal block 
consists of several 4 x 4-blocks. The number of 4 x 4-blocks is determined by the order 
of finite elements: (uord + 1)^- From Eq.(|5^ we notice that the occupation pattern in a 
4 X 4-block results from a superposition of occupation patterns of the matrices (56). 
Non-zero matrix elements outside the block diagonal band structure result from nonlocal 
terms that correspond to the finite range pairing interaction (dark shaded squares in Fig. 
3 (a)). In the Figs. 3 (c) and (d), and Figs. 3 (e) and (f) we display the occupation pat¬ 
terns for the choice of second and third order finite elements, respectively. To each pair 
of nodes {p,p') on the finite element mesh, there corresponds a 4 x 4-block {wp\H\Np>), 
which has to be multiplied by a vector Xpi. Depending on the boundary conditions (|75|), 
one or more componets of Xpi are set to zero. Boundary conditions are taken into account 
by simply eliminating columns and rows with the corresponding index from both stiffness 
matrices A and N. We illustrate this procedure in Figs. 3 (a) and (b) for the case k = —1 
(boundaries r = 0 and r = r ma x). In the Figs. 3 (c) and (d) the boundary conditions 
correspond to k = -|-1, while in the Figs. 3 (e) and (f) boundary conditions for the case 
|k| > 1 are displayed. 


4 An illustrative calculation 


In this section we present an illustrative calculation for the ground state of the spherical 
nucleus ^^^Sn. The self-consistent calculation is performed for the mean-field parameter 
set NL1||^, and the DIS |1£] parameters for the finite range pairing interaction (|3^. The 
Z = 50 protons in Sn form a closed shell. Single-particle wave functions of proton states 
are calculated as solutions of the radial Dirac equation, as described in Ref. |T^. N = 74 
neutrons in the isotope ^^^Sn partially occupy the major shell N = 50 — 82. For neutrons 
we solve the radial RHB-equations (^6|), and Klein-Gordon equations (|^ to (|40|) for the 
meson and photon fields. In the initial step of the self-consistent iteration, Woods-Saxon 
shapes are used 


S{r) = S’(0)(l -kexp (^^^——)'j 


-1 


(77) 







(78) 


V{r) = V (0) (^1 + exp (~^)) 

for the scalar a and vector oj potentials. The contribution of the p-meson to the effective 
potential is set to zero in the first iteration step. The initial Coulomb potential corresponds 
to a homogenous spherical charge distribution of radius r*. For ^^^Sn the initial values 
of the scalar and vector potentials at r = 0 are chosen 5(0) = —395MeV and F(0) = 
320 MeV, respectively, a = 0.5 fm and rg = 6.0 fm. The RHB equations are solved for 
K = ±1, ...±7. The inclusion of additional K-blocks in the self-consistent calculation did 
not change the results for ground state properties. As is illustrated in the energy diagram 
in Fig. 2, the diagonalization of the resulting eigenvalue problem is performed in a window 
of positive quasi-particle energies [0,100] MeV. 

In each step of the self-consistent iteration, the value of the chemical potential A has to be 
adjusted in such a way that the expectation value of the particle number operator equals 
the number of neutrons, i.e. 74 in this case. The correct value of A is found as the root 
of the function 


OO 

dn{X) := (79) 

ur 'n 


where N is the number of neutrons, and the second term is the trace of the density 
matrix. In the initial steps of the self-consistent iteration, A is not calculated with a 
very high precision. The precision increases with the number of iteration steps, i.e. with 
the accuracy with which the mean-field is calculated. If Aq is the calculated value of 
the chemical potential in the i-th iteration step, we define an intervall := [Aq^ — 

AX^^\ Aq^ -|- AA^*)] in which the new value Ag"*"^^ is to be found. The width of the interval 
is dehned 




AA(o) 


(80) 


where AA^^^ is an input parameter, and we take n = 2. This procedure leads to good 
convergence, and accurate values of A are obtained when the iterations approach the self- 
consistent solution. After the first few RHB iterations, the number of nested A-iteration 
steps varies between 1 and 3. 

In what follows we present results for the self-consistent solution that correspond to 
the ground state of The meson fields have been calculated with precision 10“^ MeV, 

and the accuracy for single quasi-particle energies is 10“'^ MeV. In Figs. 4a, 4b, 4c, and 
4d we display the normalized amplitudes of the quasi-particle spinors. The normalization 
for quasi-particle states is defined: 


1 = J -h 4>|^(r)4>v'(r)) (81) 

In Fig. 5 we display the baryon densities (jH) calculated with 150 linear shape functions 
in a radial box of 30 fm. The dashed curve corresponds to the density calculated in 
the first iteration step with the initial Woods-Saxon potentials, and the solid curve is 
the the self-consistent result. Since the densities are very sensitive to the details of the 
numerical approximations, we use them to choose the most effective cut-off in the stiffness 
matrices of the Hartree-Bogoliubov equations. Namely, as explained in Section 5, we only 
construct stiffness matrices with band diagonal structure. The width of the bands is an 




input parameter, i.e. the cut-off parameter of the matrices. A too small value of the cut¬ 
off parameter means that many off-diagonal matrix elements of the pairing interaction are 
neglected. A large cut-off makes the diagonalization of the algebraic eigenvalue equations 
slow. In the calculation of the baryon density in Fig. 5 we have used a cut-off value of 50. 
The resulting global stiffness matrices have a band diagonal structure of width 101. This 
value is too small, and unphysical oscillations of the density are observed. The amplitude 
of these oscillations corresponds to the values of the largest pairing matrix elements that 
have been neglected. Most of our calculations, and in particular the results that follow, 
have been performed with a cut-off parameter 60. 

In Fig. 6 we show, for each value of the K-quantum number, the contributions of 
individual quasi-particle states to the pairing field. We plot the quantities 




\ 0 (r, r') 


jiVir') 


(82) 


The contributions to the pairing field are mainly concentrated on the the surface of the 
nucleus. In Fig. 7 we also display the integrated kernel of the integral operator of the 
pairing interaction 


/A(s)(r)\ 
(r )) 


'f' max 

J dr'r''^ 
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A(99)(r, r') 
0 


A(//)(r, r') ) ( 1) 


(83) 


In order to illustrate how the single-particle density and pair matrices depend on the size 
5"max of the spherical box, in Fig. 8 we display the quantities Nnn and Pnn 

max 

NnK= J drr^{gZ\r)gl^J{r) + (84) 

0 


PnK = - J drr^{g^^J{r)g22ir) + (85) 

0 


which constitute a measure of the contribution of the n-th quasi-particle state to the 
density matrix p and to the pairing tensor k, respectively. {2j + l)Nn is the contribution 
to the particle number. In Figs. 8a and 8b we plot Nn and Pn as functions of quasi-particle 
energy for the Si/ 2 -states (k = —1). Self-consistent results are shown for three sizes of the 
radial box: rmax = 15, 30 and 40 fm. The precision of the calculations is: 10“^ MeV for 
the meson helds, and 10“^ MeV for the quasi-particle energies. The largest contributions 
to the density p come from the three quasi-particle states at 1.52 MeV, 26.18MeV and 
51.82MeV). The state at 1.52MeV, which is closest to the Fermi level A, gives the largest 
contribution to the pairing tensor. For k = — 1 the contribution to the pairing comes from 
these three peaks that represent bound states. Both quantities, iV„ and Pn-, asymptotically 
decay with increasing quasi-particle energy. Therefore, an extension of the energy window 
above 100 MeV would not change the calculated quantities. Similar results are obtained 
for other values of k. In Tables 1 and 2 we illustrate the numerical accuracy with which 
calculations can be performed. From 80 to 200 linear finite elements have been used in 
the calculations. 14 values of k are included: k = —7 to k = +7. In Table 1 quantities 
that characterize the bulk properties of ^^^Sn are displayed. For all quantities we observe 
convergence with the increase of the size of the radial box, and of the quasi-particle energy 


window. In addition, in the right column, we list the corresponding values calculated with 
a computer program that uses an expansion in oscillator basis functions . For a nucleus 
that is still far away from the drip line, the finite element discretization and the oscillator 
basis expansion should produce very similar results. For drip line nuclei we expect the 
finite element method to provide more accurate solutions, rmax and F^max denote the size 
of the radial box and of the energy window, respectively. The energies Ep and E„ are 
defined as 

^n(v) — ^\^i\'^i.n.fp)-^i.n(v)j (^ 6 ) 

i 

where the denotes the canonical energies, and the occupation probabilities. 

The pairing field Aaair,r') and the pairing tensor Kaair,r') are used to calculate the 
pairing energy Epair 


Epair — / dvv 


oo 

J drr'^'^Kaair, r')Aaa{r, r'). 


(87) 


The sum in (|87| ) runs over all values of k. In addition we display in Table 1 the total 
energy of the u-held Egig, the contribution of nonlinear terms to the u-energy Eni, the 
energy of the w-field Eomej the energy of the p-field Erho, the Coulomb energy Epho, the 
total binding energy Eb and the binding energy per nucleon These quantities are 


defined in Refs. |]^| and |£3. The mean square radii rms are dehned 


rms = 


\ 


0_ i 

oo 

/ drr2X;2|Ki|d>[, .(r)d>y,i(r) 


( 88 ) 


In Table 2 we compare results for canonical energies and occupation probabilities of single¬ 
neutron states in 


5 Program structure 


The program is coded in C++. The implementation of the relativistic mean field model 
in the Hartree approximation for spherical doubly-closed shell nuclei has been described 
in Ref. |]^]. In this section we only describe the changes that have been made in order to 


extend the program to open-shell nuclei, and include pairing correlations in the framework 
of relativistic Hartree-Bogoliubov theory. 

The main part of the program consists of seven classes: MathPar: numerical pa¬ 
rameters used in the code, PhysPar: physical parameters (masses, coupling constants, 
etc.), FinEl: finite elements. Mesh: mesh in coordinate space. Nucleon: neutrons and 
protons in the nuclear system. Meson: mesons and photon with correponding mean 
helds and the Coulomb field, and the class LinBCGOp. A detailed description of these 
classes can be found in Ref. [I3|. In the following we describe additional methods and 


parameters that are used in the present version of the program. The program allows one 
type of nucleons (protons or neutrons) to partially occupy a major shell of the nuclear 
shell model. A partially occupied major shell is referred to as an open shell. Nuclei with 
both proton and neutron shells open are generally deformed. The program is restricted 








to nuclear systems with spherical symmetry. For the type of nucleons that occupy open 
shells the relativistic Hartree Bogoliubov (RHB) equations have to be solved for a set 
of K-quantum numbers: ±1, ±2, ±3, ... The class MathPar contains the new parame¬ 
ter Kapa-max which equals the maximal absolute value of k. Other new parameters are: 
Mesh_rmax_mes for the size of the mesh on which the meson field equations are solved, Nu- 
cleon-maxrhbst[PhysPar::NumNucTypes] for the maximal number of quasi-particle states, 
LambdaO for the initial value of the chemical potential, D_LambO for the initial step in the 
A-iteration procedure, and Numb_rhb_tol for the initial tolerance in the calculation of the 
number of particles for a particular value of A. Due to nonconservation of the number of 
particles, solutions of RHB equations in general do not correspond to the correct number 
of nucleons. In the A-iteration the value of the chemical potential is adjusted in such 
a way that the expectation value of the particle number operator (trace of the density 
matrix), equals the number of nucleons for the specific nucleus. 

The set of physical parameters of the nucleus, contained in the class PhysPar, is ex¬ 
tended by the parameter RHB_nuc-number[PhysPar::NumNucTypes], which is the number 
of nucleons for which the RHB equations are solved. The prototype method set-nucleus_XX() 
can be used to define the numbers of protons and neutrons. In addition to the mean field 
parameter sets NLl, NL2, NLSH and NL3, the class PhysParcontains two parameter 
sets for the Gogny interactio: D1 |l^ and DIS |I^] . 

Essential changes have been made in the class nucleon. The method solve( Meson 
constSz _sigma, Meson constX, -omega, Meson eonstX, -rho, Meson const&L -photon ) dis¬ 
tinguishes between nucleons that occupy closed shells, for which the Dirac equation has 
to be solved, and nucleons in open shells, for which the wave functions are solutions 
of Hartree-Bogoliubov equations. A new constructor make-gog-stiff(int nvec,int ikap,int 
*iab,int *iz,double* A,double *N) for stiffness matrices of the RHB equations has been 
added to the class, make-gogstiff is used only in the first step of the A-iteration, for each 
K. The method store-gog-stiff(int nvec,int ndx,int ikap,int* iab,int* iz, double* A, dou¬ 
ble *N). stores the resulting stiffness matrices in the private fields int * *siz, int * *siab, 
double**SA, and double**SB. If the diagonalization has to be repeated just for a different 
value of the chemical potential, the reconstruction of the stiffness matrices is performed by 
read-gog-Stiff(nvec,ndx,ikap,iab,iz,A,N) and make-newJamb-gogstifffA,lambda,lamb-old). 
Matrix elements that depend on A are replaced with new values. 

houeond(kapa,nvec,n-eig,iab,iz,A,N,k.nvec2,k,n2) includes boundary conditions in the stiff¬ 
ness matrices. Rows and columns are eliminated, and the resulting matrices are written 
in condensed form. 

The largest matrix elements of the finite range pairing interaction are concentrated 
near the diagonal of the stiffness matrix. The absolute values of the pairing matrix el¬ 
ements monotonically decrease as one goes away from the diagonal, and the outermost 
are many orders of magnitude smaller than those at the diagonal. Therefore, elements of 
the stiffness matrices are only calculated in blocks within a band diagonal structure. The 
mat-eutoff(int nvec2,int n2,int* iz,int* iab,double *A,double* N, int* nvec3,int* ndfS) 
removes from the global stiffness matrix all elements which lie outside a band of chosen 
width, and performes a condensation of the matix. The eigenvalue problem with band ma¬ 
trix structure is solved by the bisection method bisec(n2,A,N,iab,iz,en-low-rhb,en-upp-rhb,tol, 
lam,XX,ndx,mdx,ndf3). Eigenvalues and eigenvectors are calculated in a window 
\low-rhh-, upp-vhh] of the energy parameter. 

The function elim-spur-rhb(mdx,XX,lam,ist,k,numb-lev,kapa) eliminates spurious solu¬ 
tions from the spectrum produced by bisee. add-dens-rhb adds contributions to the source 




terms of the boson fields. For a calculation performed with some value of the chemical 
potential, the method trace-dens(int ist,int numbJev,int kapa) calculates the trace of the 
density matrix for each k. A combination of the secant method and the bisection method 
is used by eale-newJambda-secant(lamb-iter,dn,SzlaJ,hla-r,k,dnJ,Szdn-r, lam_old) to cal¬ 
culate new A values. In the following table we list all the new member functions which 
have been included in the class nucleon. 


List of new member functions in class nucleon: 

void write_canon_energies( double** vv, double** Ec ); 

void write_solutions(int kapa,double *lam,double **XX,int mdx ); 

void write_rhb_states(int k); 

void write_canonicaLbasis(int kapa,int n,double **CG,double **CF); 
void write_dens_kap(int k); 

void write_stiff_matels(int k,int nvec,int n,int *iz,int *iab, double *A,double *B); 

void calc_canonicaLbasis(); 

void calc_NnPn(); 

void scan_delta_wave(int k); 

void scan_delta(); 

double energy_pair(); 

double energy_pair_delta(); 

double radius_ms(); 

double sum_r2(); 

double norm(); 

double fac(double n); 

void make_kapa_list(); 

double guvwave( double const* g, int kapa, int ife, int iga ) const; 
double fuvwave( double const* f, int kapa, int ife, int iga ) const; 

void make_stiff_delta_force(int lamb_iter,int nvec,int ikap, int *iab,int *iz, double *ast,double 
*bst); 

void make_gog_stiff(int nvec,int k,int *iab,int *iz, double *ast,double *bst); 

void make_newJamb_gog_stiff(double *ast,double lamb_new, double lamb_old); 

void store_gog_stiff(int nvec,int n,int ikap,int *iab,int *iz, double *A,double *B); 

void read_gog_stiff(int nvec,int n,int ikap,int *iab,int *iz, double *A,double *B); 

void boucond(int kapa,int nvec,int n,int *iab,int *iz, double *ast,double *bst,int *nvec2,int 

*n2); 

void mat_cutoff(int nvec,int n,int *iz,int *iab, double *ast,double *bst,int *nvec3,int *ndf3); 
void calc_gofac(); 

void make_s3_locst( int ife, double** s3 ); 

void make_loc_gog_st(int ikap,int ifel,int ife2,double** s6, double** s7); 
void make_loc_stiff_delta(int ikap,int ifel,double **s6, double** s7); 
double delta(int swgf,int kl,int ifel,int igal,int ife2,int iga2); 
double delta_delt(int swgf,int kl,int ifel,int igal); 

void elim_spur_rhb(int mdx,double** XX,double const* lam, int ist,int *numb_lev,int kapa 

); 

void pickup_rhb_state(int iw,int kapa,double const* X,double lam); 



void setup^tiffmats(); 

void copy_new_old(); 

void calc_densO_rhb(); 

double density_rhb(int ife,mt iga); 

void make_dens_matr(double *A,double *B,mt *iab,int *iz, int cutoff,mt *nvec,int *nn); 
void make_kk_dens_matr(int k,double *A,double *B,int *iab,int *iz, int *nvec,int *n); 
void make_dens_matr_sdiag(int k,double **A,int *n); 
void add_dens_rhb(int ist,int numbJev,int kapa); 
void large_r_dens (double r_cutoff); 

void bisec_lambda(int lambJter,double dn, double *la_l,double *la_r, double *dn_l,double 
*dn_r); 

void bisec_lambda(int lambJter,double dn,double *dnl,double *dn2); 
void calc_newJambda_BCS(); 

void calc_newjambda_secant(int lambJter,double dn, double *laJ,double *la_r, double 
*dnJ,double *dn_r, double *lam_old); 
double calc_nfe_mes(); 

double pairJens_kap(int swgf,int k,int ifel,int igal, int ife2,int iga2); 

double trace_dens(int ist,int numb Jev,int kapa); 

double deltaJambda(double dn,double sum); 

double add^um(int ist,int numb Jev,int kapa); 

double gofacO(double ll,double 12,double jl,double j2,double lam); 

double gofacl(double ll,double 12,double jl,double j2,double lam); 

double pl(int 1, double x); 

double vl(int l,int ifel,int igal,int ife2,int iga2,double mu_gog); 

double clebsh_gordon_coeff(double jl,double ml, double j2,double m2, double j3,double 
m3); 

double wig_3j(double jl, double j2, double j3, double ml, double m2, double m3); 
double wig_6j(double jl, double j2, double j3, double 11, double 12, double 13); 
double del(double jl, double j2, double j3); 

For the meson fields there is a possibility to choose the size of the finite element mesh dif¬ 
ferent from that used in the solution of the Hartree-Bogoliubov equations. In some cases 
we have found that, by taking a smaller radial box for the meson fields, better numerical 
accuracy is obtained for quasi-particle spinors that correspond to nucleon states in the 
continuum. The number of finite elements in the meson mesh is determined by the mem¬ 
ber function double Meson::calc-nfe-mes() in the constructor of the class meson. The 
size parameter of the radial box r_max-mes has been included in the class mesh. In nu- 
mutils.cc we have included a new solver for matrix diagonalization void sdiag(int n, double 
**a, double *d, double **x, double *e, int is) . It is based on the Householder algorithm, 
and is used for matrices of dimension smaller than 1000 x 1000. sdiag diagonalizes the ex¬ 
tremely ill conditioned density matrices pkk' in the member function calc_canonicaLbasis 
of class nucleon. The diagonalization is performed in the final transformation of single 
quasi-particle solutions to the canonical basis of single particle states. 



A RHB equations for nuclear systems with spher¬ 
ical symmetry 


The coordinate transformation 

T : [0, oo) (8) [0, tt) (g) [0, 27r) —> 

I—> (r sin 0 cos (/), r sin 0 sin 0, r cos 0) (89) 

defines the RHB equation (|2[) in spherical coordiantes {r,9,4>). 

Hd = -i (as dr + ai-de + 02 — d^) + V{r) + m/3 + S{r) (5 (90) 

r r sin 0 

is the Dirac hamiltonian in spherical coordinates. The standard representation is used for 
the matrices a* (i = 1, 2, 3): cii = ui (g) Uj, where ai are the Pauli matrices. For a system 
with spherical symmetry, the scalar and the vector potential read 


S{r) = gacr{r) 

y{r) = + QpT p{r) + A^{r) (91) 


Single-particle eigenfunctions of the Dirac hamiltonian Hq are at the same time eigen¬ 
vectors of the total angular momentum J, its z-component Jz, and of the operator 

k := -i(5{T.2de - (92) 

sm 9 

where Sj := I 2 <g) (i = 1, 2,3). For the eigenfunctions of H£) 

= (k = ±1,±2,±3,...), (93) 

and this relation motivates the ansatz 

0, (j)) := ignK{r),i fn-kk)'^ ^ {^K,m{9, ({>), (/>))^. (94) 

The Dirac hamiltonian takes the form 


Hd = -iiasdr + ^ 3 r ^ K) + V{r) + m {P - I 4 ) + S{r) P, (95) 


and the left hand side operator of the RHB equation (P) for a spherical symmetric system 
reads 


Hrhb = <73 (g) [—i [a^dr -I- 13 -H) -|- H(r) -|- mP -|- S{r)P] - Xa^ (g) I 4 -|- di (g) A(r). (96) 
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The ansatz for the nucleon spinor in (^) is 

/ gy'\r)^,,^rn{9A) \ 

(U\^ 

Uy g^y{r)^,,m{9,P) 

Vi/(^)(r)H_K,m(^'»/ 


(97) 


on 







where the radial functions g{r) and /(r) depend on the principal quantum number n, and 
on K. The pairing tensor is defined Kcd{r,f^) ■= and we postulate 

l^nn'kk'mm' (r,r , 9, 9 , (j), (j) — 




0 




i 0 

For a pairing interaction of finite range, the kernel Aab(r, r') (pi],) is written in the form 

A^^'{r,9,4>,r',9',cj)') = 


1 


EEE ^ )^n,h'^kK' ^mrh' ^nh^kk'mrh' (^; ^ ^ ); 


(99) 


n,n' h^k' m,m' 


and the integral operator (p!^)of the pairing field, as a function of spherical coordinates 

CX) TT 27 r 

A{r,9,(j)) = J J J dr'r'‘^d9sm9d(l)^'^V^,^,iiiiKiik{r,9,(t>,r',9',4>'). (100) 


0 0 0 


B Two-nucleon matrix elements of the Gogny in¬ 
teraction 


For isospin T=1 pairing, i.e. pairing interaction between identical nucleons, the Gogny 

( 101 ) 


force (^) can be written in the form 


V{1,2)= ^i(|ri-r2|) (Ai + BiP^^) 

i=l,2 

where the spin operator acting on the two-nucleon state is defined 


P"l(ll)S) = (2S-l)|(il)S) = { 


-1-1 for 5 = 1; 
-1 for 5 = 0; 


The radial interaction 


(ri -r2)'^ 


Vii\ri - r2|) = e 


( 102 ) 


(103) 


can be written as an infinite sum of terms, each separable in the angular coordinates of 
the two nucleons 


where 


and 


^(|ri - r2|) = E rx(n,r2)-PA(cos6»i2), 

A=o 

Att _ 

Pa(cos6'i 2) = ^^_^^ E^AM(ri)U/.(h), 


^A(n,?’2) = 


2A-M 


dcos0i2’F(|ri - r2|)PA(cos0i2) 


(104) 

(105) 

(106) 


-1 


oi 







The non-antisymmetrized, two-nucleon matrix element between non-normalized states 
reads 

VjM = {rljyi'j'\V{\v - v'\)\rllrTj')jj^, (107) 

where J and M are the total angular momentum of the two-nucleon state, and its 2 - 
projection, respectively. Since the interaction contains the spin operator, it is convenient 
to calculate the matrix element in the LS-coupling scheme. The transformation between 
the jj and LS-coupling schemes for the two nucleon state 

\{l,l'),jj-JM)=Y,jj'LS\l r l\\{11 ')SL-JM) (108) 

iS [j f J) 


where L is the total orbital angular momentum, and S the total spin of the two-nucleon 
state. We use the short-hand notation j := y/2j + 1. The matrix element can then be 
written 


VjM = EEl j'LsTfLS{rr', {ll')SLJM\V{\r - r'\)\ff', {ll')SLJM) 

LS L,S 



1 

2 

C \ ( 1 

1 1 2 

1 

2 


1 ' 


i 


1 j 

j’ 

J) 

1 j 

f 


Using the Slater expansion (1^), the LS-coupling matrix element takes the form 

A frr 

= E ^^^A(r,rOE((/Z')5^^M|y,;(f)n^(r')(^ + BP")|(/T')5LJM), (110) 


where the radial factors V\{r,r') are defined in (|106| ), and the angular part is calculated 
j:{{inSLJM\Y^r)Yxf,i?){A + Bp-)\{ll')SLJM) = 


6ssSrL{A + B{2S-l)){-lp+^+^ ® 


I A I 
0 0 0 


I' A [' 

0 0 0 


I V L 

I' I A 


( 111 ) 


For the pairing interaction we only consider matrix elements between two-nucleon states 
with total angular momentum J = 0. In this case the 9j-coefficients in ( |109| ) reduce to 
6j-coefficients, and the two-nucleon matrix element reads 

U7=o = 6jj,5y.,5u,5ii,{2l + l){ 2 l + l)fj' 

^ (_i)i+i+i+5(25 + i)|f [ <51 n I 

5 = 0,1 I 2 2 3 ) V 2 2 3 ) 

EU(r,/)(.4 + B(25-!))(( J ()'{( ; 1]. (112) 

The allowed integer values for A in the sum are 

\l-l\<\<l + l (113) 


The matrix elements essentially consist of two terms: 5 = 0 and 5 = 1. Introducing the 
explicit expressions for the 3j and 6j-coefficients, the two terms read 


Vj=o s=o = {A- B)5jy5y-,5iv5ji,jj Y. 

Aeven 


{2g-2l)\{2g-2X)\{2g-2l)\ 

2 ( 29 + 1 )! 


A 


2 


{g-l)\{g-\)\{g-l)\ 


(114) 










and 


Vj=o s=i = + 


(|+d^+i)—i(i+i))(|+d^'+i)-i(i+i)) 


2V- I i(z+i)i-(z+i) ' ( l)^'Kx(r,r') 


|^/(/ + 1) + l(J + 1) + A(A + 1) 
where 2g = I + X + 1. 


(2g-2l)l(2g-2X)H2g-2l)l 

2(2s+l)! 


Aeven 


a'- 


n2 


{g-iy.{g-xy.{g-l)\ 


(115) 
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Figure Captions 

Fig. 1 Vector space homomorphisms between the general single-particle basis, the ba¬ 
sis of quasi-particle states, and the canonical basis of single-particle states. The 
diagrams illustrate the transformations (^) to &■ 

Fig. 2 Relativistic Hartree-Bogoliubov model for a finite nucleus. Single-particle eigen- 
spectrum of a Dirac hamiltonian (left and center), and single quasi-particle eigen- 
spectrum of the Hartree-Bogoliubov equations (right). 

Fig. 3a Occupation pattern of the global stiffness matrix A (^) for n = —1 and linear 
shape functions. The dark grey squares indicate matrix elements of the finite range 
Gogny interaction. Due to the nonlocal character of the interaction, the matrix 
elements are distributed over the whole matrix. The light grey squares correspond 
to matrix elements which result from the local operators in (|3^). They form a block 
diagonal band structure. Boundary conditions are taken into account by eliminating 
the corresponding rows and columns. The matrix is symmetric and we display the 
global index of nonzero matrix elements, as used by the bisection method in the 
solution of the eigenvalue problem. 

Fig. 3b Occupation pattern of the overlap matrix N (^), for k = — 1 and linear shape 
functions. 


Fig. 

Fig. 

Fig. 

Fig. 

Fig. 


Fig. 

Fig. 

Fig. 

Fig. 


Fig. 


Fig. 


Fig. 


3c Same as in Fig. 3a, but for k = +1 and quadratic shape functions. 


3d Same as in Fig. 3b, but for k = -|-1 and quadratic shape functions. 
3e Same as in Fig. 3a, but for |k| > 1 and cubic shape functions. 

3f Same as in Fig. 3b, but for |k| > 1 and cubic shape functions. 


4a Self-consistent normalized quasi-particle wave functions. Radial amplitudes of 
U-components (upper left), and V-components (upper right), for the lsi /2 (solid), 
2 si/ 2 (dashed), and 3si/2 (doted) bound states in ^^“^Sn. The radial amplitudes of 
continuum si/ 2 -states are shown in the lower part of the figure: U-components (left) 
and V-components (right). The self-consistent solution is calculated with 100 linear 
shape functions on a uniform radial mesh [0, 30 fm]. 

4b Same as in figure 4a, but for pi/ 2 -waves. 

4c Same as in figure 4a, but for p 3 / 2 -waves. 

4d Same as in figure 4a, but for /i 3 i/ 2 -waves. 

5 Baryon density of ^^^Sn after the first iteration step (dashed), and for the self- 
consistent solution (solid). The results correspond to a calculation with 100 linear 
shape functions on a uniform radial mesh [0, 30 fm]. 


6 K = ±1,±2,±3 contributions to the pairing field. For each value of k we display 
the largest contributions from individual states. The quantities that we plot are 
defined in (^). 

7 The components A*'^^(r) and A(-^^(r) of the self-consistent pairing field (^), cal¬ 
culated on a radial mesh of 15fm (solid) and 30fm (dot dashed). The pairing field 
is concentrated on the surface of the nucleus. 


8 Logarithmic plots of iV„ (|3) and Pn (|85|)), as functions of quasi-particle energy for 
the Si/ 2 -states. Results are displayed for three sizes of the radial box: rmax = 15 fm 
(dashed), rmax = 30fm (dot dashed) and r^ax = 40fm (solid). 



test run 

1 

FEM 

2 

FEM 

3 

FEM 

4 

FEM 

5 

EEM 

6 

EEM 

7 

Osc. basis 

^max 

[fm] 

15.00 

20.00 

30.00 

40.00 

20.00 

20.00 


Ejiiax 

[MeV] 

100.00 

100.00 

100.00 

100.00 

70.00 

110.00 


Ep 

[MeV] 

-1320.74 

-1320.68 

-1320.79 

-1320.79 

-1320.51 

-1320.80 


En 

[MeV] 

-1719.42 

-1646.68 

-1719.56 

-1719.56 

-1719.17 

-1719.54 


Epair 

[MeV] 

-19.948 

-19.9494 

-19.9510 

-19.9342 

-19.5666 

-19.9901 

-19.68 

Esig 

[MeV] 

17126.70 

17126.20 

17126.20 

17126.2 

17125.30 

17126.40 

17432.33 

Enl 

[MeV] 

-324.004 

-323.988 

-323.983 

-323.983 

-324.019 

-323.977 

-327.910 

Eome 

[MeV] 

-14699.40 

-14699.27 

-14699.10 

-14699.00 

-14698.10 

-14699.20 

-14674.06 

Erho 

[MeV] 

-55.460 

-55.455 

-55.448 

-55.448 

-55.449 

-55.449 

-55.529 

Epho 

[MeV] 

-366.072 

-366.072 

-366.073 

-366.073 

-366.066 

-366.075 

-365.71 

Etot 

[MeV] 

-1060.58 

-1060.69 

-1060.83 

-1060.82 

-1055.37 

-1060.86 

1051.96 

EfeM 

[MeV] 

-8.553 

-8.554 

-8.555 

-8.555 

-8.550 

-8.555 

-8.48 

A 

[MeV] 

-6.65727 

-6.65718 

-6.65706 

-6.65668 

-6.65719 

-6.65729 

-6.65786 

rmSneutron 

[fm] 

4.90853 

4.90869 

4.90874 

4.90874 

4.90825 

4.90871 

4.91073 

rmSppoton 

[fm] 

4.59741 

4.59745 

4.59751 

4.59741 

4.59738 

4.59740 

4.60187 

rmSnuQieus 

[fm] 

4.78555 

4.78561 

4.78564 

4.78564 

4.78515 

4.78562 

4.78859 


Table 1 Bulk properties of Calculated quantities are compared for various sizes of the radial 

mesh, and with results of an expansion in the oscillator basis |^. 
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Tmax [ fm ] 

15.00 

20.00 

Emax [MeV] 

100.00 

100.00 

Els,/, [MeV] 

- 51.834 

- 51.834 

E 2 s,/, [MeV] 

- 39.294 

- 39.298 

E 3 s,/, [MeV] 

- 8.227 

- 8.231 

Els,/, [MeV] 

12.988 

12.997 

Eip 3 /, [MeV] 

- 43.078 

- 43.076 

E 2 p 3 /, [MeV] 

- 27.023 

- 27.029 

E3p 3/, [MeV] 

2.977 

2.969 

Eip,/, [MeV] 

- 42.976 

- 42.998 

E 2 p,/, [MeV] 

- 24.499 

- 24.489 

Esp,/, [MeV] 

3.863 

3.850 

Eid 5 /. [MeV] 

- 38.359 

- 38.358 

E 2 d 3 /, [MeV] 

- 10.678 

- 10.683 

E3d3/, [MeV] 

13.245 

13.265 

Eid 3 /, [MeV] 

- 35.690 

- 35.691 

E 2 d 3 /, [MeV] 

- 8.489 

- 8.494 

E3d3/, [MeV] 

13.851 

13.887 

Eig,/, [MeV] 

- 16.971 

- 16.970 

E 2 g,/, [MeV] 

11.400 

11.512 

Eih,,/, [MeV] 

- 6.349 

- 6.348 

E 2 h„/, [MeV] 

22.127 

22.116 

^lSl/2 

0.999995 

0.999995 

^is,/. 

0.999598 

0.999596 

V? 

3s,/, 

0.984564 

0.984573 

^4si/2 

0.000048 

0.000048 

V? 

lP3/2 

0.999900 

0.999900 

^2p3/2 

0.998427 

0.998419 

^3p3/, 

0.000499 

0.000503 

"^lpi/2 

0.999941 

0.999940 

"^ip,/. 

0.999025 

0.999019 

v2 

3Pl/2 

0.000218 

0.000219 

"^W5/2 

0.999205 

0.999200 


0.987226 

0.987186 

'^ids/. 

0.000135 

0.000135 

< 

0.999384 

0.999381 

< 3/2 

0.966344 

0.966348 

< 3/2 

0.000095 

0.000095 

< 9/2 

0.992754 

0.992711 

< 9/2 

0.000410 

0.000415 


3 4 5 6 


30.00 

40.00 

20.00 

20.00 

100.00 

100.00 

70.00 

110.00 

- 51.834 

- 51.836 

- 51.828 

- 51.832 

- 39.303 

- 39.300 

- 39.290 

- 39.305 

- 8.235 

- 8.235 

- 8.221 

- 8.235 

13.050 

12.902 

12.981 

13.085 

- 43.075 

- 43.077 

- 43.069 

- 43.067 

- 27.035 

- 27.032 

- 27.018 

- 27.043 

2.961 

2.961 

2.983 

2.965 

- 43.006 

- 43.008 

- 42.923 

- 43.001 

- 24.481 

- 24.479 

- 24.498 

- 24.486 

3.843 

3.842 

3.894 

3.849 

- 38.357 

- 38.357 

- 38.358 

- 38.359 

- 10.686 

- 10.686 

- 10.671 

- 10.685 

13.271 

13.218 

13.203 

13.283 

- 35.692 

- 35.692 

- 35.681 

- 35.693 

- 8.496 

- 8.497 

- 8.471 

- 8.496 

14.051 

13.801 

13.766 

14.096 

- 16.969 

- 16.969 

- 16.998 

- 16.969 

11.557 

11.551 

11.289 

11.466 

- 6.347 

- 6.347 

- 6.355 

- 6.347 

22.107 

22.110 

22.138 

22.212 

0.999995 

0.999995 

0.999995 

0.999995 

0.999594 

0.999594 

0.999611 

0.999593 

0.984585 

0.984603 

0.984593 

0.984558 

0.000049 

0.000049 

0.000046 

0.000049 

0.999899 

0.999899 

0.999835 

0.999899 

0.998413 

0.998414 

0.998483 

0.998410 

0.000505 

0.000504 

0.000413 

0.000506 

0.999939 

0.999939 

0.999964 

0.999939 

0.999015 

0.999016 

0.999082 

0.999013 

0.000220 

0.000220 

0.000203 

0.000221 

0.999198 

0.999198 

0.999256 

0.999195 

0.987160 

0.987166 

0.987347 

0.987134 

0.000136 

0.000135 

0.000128 

0.000137 

0.999379 

0.999378 

0.999421 

0.999376 

0.966337 

0.966350 

0.966473 

0.966276 

0.000096 

0.000096 

0.000077 

0.000097 

0.992680 

0.992687 

0.992934 

0.992660 

0.000417 

0.000417 

0.000352 

0.000419 














